####
# Author: M. Kenwick, S. Lee, B. Kolcak
# Purpose: Perform core regression analysis of Congressional speech
# Date: May 19, 2025 
####

################################################################################
# 0 - Load Data and Packages
################################################################################
rm(list=ls())
library(lme4)
library(dplyr)
library(ggplot2)
library(stargazer)
library(readxl)

setwd('~/Dropbox/cmr_cong/replication')
data <- read.csv('master.csv')


################################################################################
# 1 - Generate and Transform Data/Variables
################################################################################
# Dropping pre-9/11 speeches
data$date <- as.Date(data$date,format="%Y-%m-%d")
data <- data[data$date >= as.Date("2001-09-11",format="%Y-%m-%d"),]

# Generate date counters
data$days <- as.numeric(data$date - as.Date("2001-09-11",format="%Y-%m-%d"))
data$months <- data$days/30

# speakerid concatenates both speakerid and session number. This splits the two
data$speaker_session_id <- data$speakerid
data$speakerid <- as.numeric(substr(as.character(data$speaker_session_id),4,9))
data$session <- as.numeric(substr(as.character(data$speaker_session_id),1,3))

# Merging in birth year and DW Nominate scores, first by integrating ICPSR IDs
key <- read.csv('members.csv') #for ICPSR IDs
data <- merge(data, key, by.x="speaker_session_id", by.y="speakerid", 
              all.x=T, all.y=F)

#Pull in Voteview data
voteview <- read.csv('voteview.csv')
voteview <- voteview %>% 
  dplyr::select(icpsr, congress, nominate_dim1, 
                nominate_dim2, born, party_code) %>% 
  unique() 
data <- merge(data, voteview, 
              by.x=c('session','icpsr'),
              by.y=c('congress','icpsr'),
              all.x=T, all.y=F)
data <- dplyr::rename(data, dwnom1 = nominate_dim1, dwnom2 = nominate_dim2)
data$age <- round(data$year - data$born,1)
data$party_code <- NULL

# Transformations 
data$tone <- (data$positive - data$negative) / data$keyword_count 
data$rep <- ifelse(data$party=="R",1,0)
data$dem <- ifelse(data$party=="D",1,0)
data$mil <- ifelse(data$military_collapse>0, 1,0)
data$mil_rep <- data$rep * data$mil
data$word_count2 <- data$word_count^2
data$ln_word_count <- log(data$word_count)
data$woman <- ifelse(data$gender=="F",1,0)
data$dw_mil <- data$mil * data$dwnom1
data$mil_cas <- data$mil * data$geo_fatal
data$mil_cas_cumu <- data$mil * data$geo_fatal_cumulat
data$mil_nocombat <- ifelse(data$mil==1 & data$combat==0, 1, 0)
data$dw_nocomb <- data$dwnom1 * data$mil_nocombat
data$dw_comb <- data$dwnom1 * data$combat
data$obama <- ifelse(data$date >= as.Date("2009-01-20",format="%Y-%m-%d"),1,0)
data$pip <- ifelse((data$obama==1 & data$dem==1)| (data$obama==0 & data$rep==1), 1, 0)
data$white <- ifelse(data$race=="white",1,0)
data$senate <- ifelse(data$chamber=="S",1,0)
data$months_rounded <- round(data$months)
data$uncertainty1 <- (data$tidy_unc1) / data$keyword_count 
data$uncertainty2 <- (data$tidy_unc2) / data$keyword_count 


# Incorporate vote choice on authorizing force in Iraq
house <- read_xls('votechoice/iraq_auth_house.xls')
senate <- read_xls('votechoice/iraq_auth_senate.xls')
auth <- rbind(house,senate)
auth$auth_yea <- ifelse(auth$V1==1,1,0)
auth$auth_nay <- ifelse(auth$V1==6,1,0)
auth <- dplyr::select(auth,icpsr,auth_yea,auth_nay)
data <- merge(data, auth, by="icpsr",all.x=T,all.y=F)
data$auth_yea[is.na(data$auth_yea)] <- 0
data$auth_nay[is.na(data$auth_nay)] <- 0
data$auth_abs <- ifelse(data$auth_yea==0 & data$auth_nay==0,1,0)

data$age_range <- NA
data$age_range[data$age < 40] <- "Under 40"
data$age_range[data$age > 39 & data$age < 60] <- "40s and 50s"
data$age_range[data$age  > 59 & data$age <80] <- "60s and 70s"
data$age_range[data$age  > 79] <- "More than 80"

#Logged national, cumulative fatalities
data$log_nat_fatal_cumulat <- log(data$nat_fatal_cumulat + 1)




################################################################################
# 2 - Descriptive Statistics (Appendix)
################################################################################
#descriptive data object
descr <- data %>% dplyr::select(tone, geo_fatal, geo_fatal_cumulat,  mil, 
                                combat, dwnom1, dwnom2, woman, age, pip, obama, 
                                senate, word_count, word_count2, months)

#Table A13 (basic descriptives)
stargazer(descr, type = 'text') 


# Table A4: War Speeches by Year (Table A4)
descr_year <- as.data.frame(table(data$year))
descr_year <- descr_year %>% mutate(Perc = Freq / 36456)
descr_year

# Table A5: War Speeches by Chamber
descr_chamber <- as.data.frame(table(data$chamber))
descr_chamber <- descr_chamber %>% mutate(Perc = Freq / 36456)
descr_chamber

# Table A6: War Speeches by Party
descr_party <- as.data.frame(table(data$party))
descr_party <- descr_party %>% mutate(Perc = Freq / 36456)
descr_party

# Table A7: War Speeces by Military Experience
descr_mil <- as.data.frame(table(data$mil))
descr_mil <- descr_mil %>% mutate(Perc = Freq / 36456)
descr_mil

# Table A8: War Speeces by Combat Experience
data$milcomb <- data$mil
data$milcomb[data$combat==1] <- 2
descr_mil2 <- as.data.frame(table(data$milcomb))
descr_mil2 <- descr_mil2 %>% mutate(Perc = Freq / 36456)
descr_mil2

# Table A9 : War Speeches by Age
descr_age <- as.data.frame(table(data$age_range))
descr_age <- descr_age %>% mutate(Perc = Freq / 36456)
descr_age 

#Individual-level data
data_ind <- data %>% dplyr::select(icpsr, state, mil, party)
data_ind <- data_ind %>% distinct()

# Table A11: Legislators Delivering War Speeches by Party
table(data_ind$party)

# Table A12: Legislators Delivering War Speeches by Military ex.
table(data_ind$mil)

#Table A10: Legislators Delivering War Speeches by Year
data_ind_year <- data %>% dplyr::select(icpsr,year, chamber) %>% distinct()
table(data_ind_year$year, data_ind_year$chamber)

#### Military Experience Descriptives in Appendix A
meta <- read.delim("fullmeta_vet.txt", header = TRUE, sep = "|", quote = "")
meta <- dplyr::filter(meta, cong>=107 & cong<=113)
meta$milcat<-"Missing"
meta$milcat[meta$military_collapse==0]<-"Non-Veteran"
meta$milcat[meta$military_collapse==1]<-"Veteran"
meta$milcat[meta$military_collapse==2]<-"Veteran"
meta$combatcat<-"No Combat Experience"
meta$combatcat[meta$combat==1]<-"Combat Experience"
meta$combatcat <- as.factor(meta$combatcat)
milcat_dist <- meta %>% group_by(cong) %>% count(milcat, combatcat)
milcat_dist$milcombat_cat <- ifelse(milcat_dist$milcat == "Non-Veteran" & 
  milcat_dist$combatcat == "No Combat Experience", "Non-Veteran",
  ifelse(milcat_dist$milcat == "Veteran" & 
  milcat_dist$combatcat == "No Combat Experience", "Non-Combat Veteran",
  ifelse(milcat_dist$milcat == "Veteran" & 
  milcat_dist$combatcat == "Combat Experience", "Combat Veteran", NA)))

milcat_dist$milcombat_cat <- as.factor(milcat_dist$milcombat_cat)
levels(milcat_dist$milcombat_cat)
milcat_dist$milcombat_cat <- factor(milcat_dist$milcombat_cat ,
                                    levels=c('Non-Veteran', 
                                    'Non-Combat Veteran', 'Combat Veteran' ))
####Figure A1
ggplot(milcat_dist) + 
  geom_bar(aes(fill = milcombat_cat, y = n, x = as.factor(cong)),
           position="fill", stat="identity",color="black",size=.4) + 
  xlab("Session") + ylab("Percent of Serving Members") + 
  scale_fill_manual(values=c("#762a83", "#d9ef8b", "olivedrab4")) +
  scale_x_discrete(breaks=c("107","108","109","110","111","112", "113", "114"),
                   labels=c("107","108","109","110","111","112", "113", "114")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank()) + labs(fill = "")

#Branch
table(meta$cong, meta$branch, useNA = "always")

branc_dat <- setNames(data.frame(matrix(ncol = 3, nrow = 49)), c("cong", "branchcat", "n"))

# army 1
branc_dat[1,] <- c(107, "army", 25)
branc_dat[2,] <- c(108, "army", 28)
branc_dat[3,] <- c(109, "army", 26)
branc_dat[4,] <- c(110, "army", 25)
branc_dat[5,] <- c(111, "army", 31)
branc_dat[6,] <- c(112, "army", 28)
branc_dat[7,] <- c(113, "army", 23+3)
# marine corps 2
branc_dat[8,] <- c(107, "marine corps", 3)
branc_dat[9,] <- c(108, "marine corps", 5)
branc_dat[10,] <- c(109, "marine corps", 3)
branc_dat[11,] <- c(110, "marine corps", 3+1)
branc_dat[12,] <- c(111, "marine corps", 5+1)
branc_dat[13,] <- c(112, "marine corps", 9+1+1)
branc_dat[14,] <- c(113, "marine corps", 8+1)
# navy 3
branc_dat[15,] <- c(107, "navy", 2)
branc_dat[16,] <- c(108, "navy", 3)
branc_dat[17,] <- c(109, "navy", 3)
branc_dat[18,] <- c(110, "navy", 4+1)
branc_dat[19,] <- c(111, "navy", 9+1)
branc_dat[20,] <- c(112, "navy", 8+1+1)
branc_dat[21,] <- c(113, "navy", 8+1+1+1)
# air force 4
branc_dat[22,] <- c(107, "air force", 10)
branc_dat[23,] <- c(108, "air force", 11)
branc_dat[24,] <- c(109, "air force", 11)
branc_dat[25,] <- c(110, "air force", 11)
branc_dat[26,] <- c(111, "air force", 10)
branc_dat[27,] <- c(112, "air force", 14+1)
branc_dat[28,] <- c(113, "air force", 13+1+1)
# coast guard 5
branc_dat[29,] <- c(107, "coast guard", 1)
branc_dat[30,] <- c(108, "coast guard", 1)
branc_dat[31,] <- c(109, "coast guard", 1)
branc_dat[32,] <- c(110, "coast guard", 1)
branc_dat[33,] <- c(111, "coast guard", 1)
branc_dat[34,] <- c(112, "coast guard", 0)
branc_dat[35,] <- c(113, "coast guard", 0)
# national guard 6
branc_dat[36,] <- c(107, "national guard", 6)
branc_dat[37,] <- c(108, "national guard", 6)
branc_dat[38,] <- c(109, "national guard", 6)
branc_dat[39,] <- c(110, "national guard", 4)
branc_dat[40,] <- c(111, "national guard", 3)
branc_dat[41,] <- c(112, "national guard", 4)
branc_dat[42,] <- c(113, "national guard", 5+3+1)
branc_dat$branchcat <- as.factor(branc_dat$branchcat)
levels(branc_dat$branchcat)
branc_dat$branchcat  <- factor(branc_dat$branchcat ,
                               levels=c('army', 'air force','national guard',
                                'marine corps','navy','coast guard'))
branc_dat$n <- as.integer(branc_dat$n)
branc_dat$cong <- as.integer(branc_dat$cong)
branc_dat_sub <- branc_dat %>% filter(branchcat!= "no service") 

####Figure A2
ggplot(branc_dat_sub) + 
  geom_bar(aes(fill = branchcat, y = n, x = as.factor(cong)),position="fill", 
           stat="identity",color="black",size=.4) + 
  xlab("Session") + ylab("Percent of Serving Members") + 
  scale_fill_manual(values=c("#d9ef8b","#8ba7ef","#a7ef8b", "#efd38b","#8bd9ef","#8fb118")) +
  scale_x_discrete(breaks=c("107","108","109","110","111","112", "113", "114"),
                   labels=c("107","108","109","110","111","112", "113", "114")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank()) + labs(fill = "")




################################################################################
# 3 - Descriptive Visualizations of Tone
################################################################################
########
# Figure 1: Congressional Sentiment in Wartime Speeches 
########
# Prep data
data2 <- data %>% 
  group_by(month = lubridate::floor_date(date, "month")) %>%
  dplyr::summarize(avg_tone = mean(tone, na.rm=TRUE), n=n()) 

pdf("plots/fig1.pdf",height=6.1,width=6.6)
par(mar=c(4,4,3,3),  oma = c(0,0,0,0),mfrow=c(1,1))
plot(data2$month,data2$avg_tone,xlab="Date",
     ylab="Average Tone of Reference to Iraq/Afg.",col="white",
     ylim=c(min(data2$avg_tone),max(data2$avg_tone)))
abline(v=as.Date("2003-03-20"),col="gray30",lty=2)
text(as.Date("2003-03-20"),1.27,"Invasion of Iraq",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2004-04-27"),col="gray30",lty=2)
text(as.Date("2004-04-27"),1.27,"Abu Ghraib Report",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2009-01-20"),col="gray30",lty=2)
text(as.Date("2009-01-20"),1.27,"Obama Inauguration",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2007-02-16"),col="gray30",lty=2)
text(as.Date("2007-02-16"),1.27,"Troop Surge Debate",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2013-12-30"),col="gray30",lty=2)
text(as.Date("2013-12-30"),1.27,"ISIS Captures Fallujah",pos=2,srt=90,col="gray30",cex=.75)
points(data2$month,data2$avg_tone,
       bg=alpha("steelblue",.4), col=alpha("steelblue",1),lwd=1,cex=data2$n/200, pch=21 )
lines(lowess(data2$avg_tone~data2$month,f=.2),col="orangered",lwd=3)
dev.off()

########
# Figure A3 Congressional Sentiment by Party
########
data_rep <- data %>% filter(rep == 1) %>%
  group_by(month = lubridate::floor_date(date, "month")) %>%
  dplyr::summarize(avg_tone = mean(tone, na.rm=TRUE), n=n())
data_dem <- data %>% filter(dem == 1) %>%
  group_by(month = lubridate::floor_date(date, "month")) %>%
  dplyr::summarize(avg_tone = mean(tone, na.rm=TRUE), n=n())

pdf("plots/figA3.pdf",height=6.1,width=6.6)
par(mar=c(4,4,3,3),  oma = c(0,0,0,0),mfrow=c(1,1))
plot(data_rep$month,data_rep$avg_tone,xlab="Date",
     ylab="Average Tone of Reference to Iraq/Afg.",col="white",
     ylim=c(min(data_rep$avg_tone),max(data_dem$avg_tone)))
abline(v=as.Date("2003-03-20"),col="gray30",lty=2)
text(as.Date("2003-03-20"),1.5,"Invasion of Iraq",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2004-04-27"),col="gray30",lty=2)
text(as.Date("2004-04-27"),1.5,"Abu Ghraib Report",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2009-01-20"),col="gray30",lty=2)
text(as.Date("2009-01-20"),1.5,"Obama Inauguration",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2007-02-16"),col="gray30",lty=2)
text(as.Date("2007-02-16"),1.5,"Troop Surge Debate",pos=2,srt=90,col="gray30",cex=.75)
abline(v=as.Date("2013-12-30"),col="gray30",lty=2)
text(as.Date("2013-12-30"),1.5,"ISIS Captures Fallujah",pos=2,srt=90,col="gray30",cex=.75)
points(data_rep$month,data_rep$avg_tone,
       bg=alpha("firebrick2",.4), col=alpha("firebrick2",1),lwd=1,cex=data_rep$n/200, pch=21 )
lines(lowess(data_rep$avg_tone~data_rep$month,f=.1),col="firebrick4",lwd=3)
points(data_dem$month,data_dem$avg_tone,
       bg=alpha("dodgerblue",.4), col=alpha("dodgerblue",1),lwd=1,cex=data_dem$n/200, pch=21 )
lines(lowess(data_dem$avg_tone~data_dem$month,f=.1),col="dodgerblue3",lwd=3)
dev.off()

########
# Figure 3: Average Speaker Tone in Wartime References by Ideology and Veteran Status
########
# Weighting the mean
data2 <- data %>% 
  group_by(speakerid) %>% 
  mutate(tone = mean(tone,na.rm=T), 
         mil = max(mil,na.rm=T),
         cas = sum(geo_fatal,na.rm=T),
         dwnom1 = mean(dwnom1, na.rm=T),
         rep = max(rep, na.rm=T),
         n = n()) %>% 
  dplyr::select(speakerid, state, tone, rep, mil, cas, dwnom1, n) %>% 
  unique()
data2$tone_mean <- mean(data$tone,na.rm=T)
data2$tone_var <- var(data$tone,na.rm=T)
data2$tone_speaker_var <- var(data2$tone,na.rm=T)
data2$tone_weighted <- ((data2$n/data2$tone_var)*(data2$tone) + (1/data2$tone_speaker_var)*(data2$tone_mean)) / 
  ((data2$n/data2$tone_var) +  (1/data2$tone_speaker_var))

pdf("plots/fig3.pdf",height=7,width=7.2)
par(mar=c(4,4,3,3),  oma = c(0,0,0,0),mfrow=c(1,1))
plot(NULL,# create empty plot
     xlim = c(-.75,1), 
     ylim = c(-1.75,1.75), 
     axes = F, xlab = NA, ylab = NA)  
grid()
box()
rep_mean <- mean(data2$tone_weighted[data2$rep==1],na.rm=T)
dem_mean <- mean(data2$tone_weighted[data2$rep==0],na.rm=T)

points(data2$dwnom1[data2$mil==0 & data2$rep==1],
       data2$tone_weighted[data2$mil==0 & data2$rep==1],
       bg=alpha("#a50026",.4), col=alpha("#a50026",.4),lwd=1,cex=.6+(data2$n[data2$mil==0 & data2$rep==1]/100), pch=21 )
points(data2$dwnom1[data2$mil==0 & data2$rep==0],
       data2$tone_weighted[data2$mil==0 & data2$rep==0],
       bg=alpha("#313695",.4), col=alpha("#313695",.4),lwd=1,cex=.6+(data2$n[data2$mil==0 & data2$rep==0]/100), pch=21 )
points(data2$dwnom1[data2$mil==1],data2$tone_weighted[data2$mil==1],bg=alpha("olivedrab3",.8), col=alpha("olivedrab4",1),lwd=1.4,cex=.6+(data2$n[data2$mil==1]/100), pch=21 )
axis(1)
axis(2)
mtext(side=1,"Speaker DW-Nominate Score",line=2)
mtext(side=2,"Average Tone of Reference to Iraq or Afghanistan",line=2.5, cex=1.1)
box()
legend(0.63,-1.32, legend = c("1 Reference","50 References","250 References"), pch = 21, 
       pt.cex = c(.61,1.1,3.1), 
       x.intersp = 1.4, y.intersp = 1.4, cex=.7,
       bty = "n", col = "black", pt.bg=alpha("gray80", 0.5))
dev.off()

########
# Figure 2: Top and Bottom 10 Average Tone of Speakers
########
# Taking Mean among speakers w/ >25 speeches
data2 <- data %>% 
  group_by(speakerid) %>% 
  mutate(tone = mean(tone,na.rm=T), 
         mil = max(mil,na.rm=T),
         cas = sum(geo_fatal,na.rm=T),
         dwnom1 = mean(dwnom1, na.rm=T),
         rep = max(rep, na.rm=T),
         n = n()) %>% 
  filter(n > 25) %>% 
  dplyr::select(speaker, speakerid, state, tone, rep, dwnom1, n, auth_yea, auth_nay) %>% 
  unique()

# Weighting the mean 
data2$tone_mean <- mean(data$tone,na.rm=T)
data2$tone_var <- var(data$tone,na.rm=T)
data2$tone_speaker_var <- var(data2$tone,na.rm=T)
data2$tone_weighted <- ((data2$n/data2$tone_var)*(data2$tone) + (1/data2$tone_speaker_var)*(data2$tone_mean)) / 
  ((data2$n/data2$tone_var) +  (1/data2$tone_speaker_var))
data2 <- arrange(data2, tone_weighted)

#Removing duplicates because of name misspelling 
data3 <- data2 %>% 
  dplyr::select(speakerid, tone_weighted, n) %>% 
  unique()
tail(data2, 20) # Lists names
head(data2, 20) # Lists names 

neg_names <- c(  "John Lewis (D-GA)",
                 "Cynthia McKinney (D-GA)",
                 "Rand Paul (R-KY)",
                 "Ron Paul (R-TX)",
                 "Bernie Sanders (I-VT)",
                 "Maurice Hinchey (D-NY)",
                 "Jimmy Duncan (R-TN)",
                 "Byron Dorgan (D-ND)",
                 "Jim McDermott (D-WA)",
                 "Ted Poe (R-TX)")

pos_names <- c( "Millender-McDonald (D-CA)",
                "Johnny Isakson (R-GA)",
                "Mark Pryor (D-AR)",
                "Henry Hyde (R-IL)",
                "Richard Lugar (R-IN)",
                "Michael Bennet (D-CO)",
                "Tom Cole (R-OK)",
                "John Warner (R-VA)", 
                "Nita Lowey (D-NY)",
                "Robin Hayes (R-NC)")

pos <- tail(data3$tone_weighted, 10)
neg <- head(data3$tone_weighted, 10)

pdf("plots/fig2.pdf",height=7,width=7.5)
par(mar=c(5,13,3,3),  oma = c(0,0,0,0),mfrow=c(1,1))
barplot(c(neg,pos),horiz=T,names=c(neg_names,pos_names),las=1,
        col=c("#3288bd","#3288bd","#d53e4f","#d53e4f","#D3D3D3",
              "#3288bd","#d53e4f","#3288bd","#3288bd","#d53e4f",
              "#3288bd","#d53e4f","#3288bd","#d53e4f","#d53e4f",
              "#3288bd","#d53e4f","#d53e4f","#3288bd","#d53e4f"),
        xlab="Average Tone (Weighted)",xlim=c(-1.25,1.25))

legend("bottomright", legend = c("Republican", "Democrat", "Independent"),
       fill = c("#d53e4f", "#3288bd", "#D3D3D3"), bty = "n")
dev.off()

################################################################################
# 4 - Regression Analysis
################################################################################
########
# Table 2: Veteran Status and Semantic Tone When Referencing Iraq or Afg.
########
m1 <- lmer(tone ~ mil + (1 | speakerid),  data=data)
m2 <- lmer(tone ~ mil + rep + (1 | speakerid),  data=data)
m3 <- lmer(tone ~ mil + rep + woman + age + 
                pip + obama + senate +
               word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),  data = data)
m4 <- lmer(tone ~ mil + rep + mil_rep + pip + obama + senate + woman + age + 
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),  data = data)
stargazer(m1, m2, m3, m4,
          keep = c('mil', 'rep', 'mil_rep', 'pip', 'obama', 'senate', 'woman',
                   'age', 'geo_fatal_cumulat',
                   'word_count', 'word_count2'), type = 'text')

########
# Table 3 (Table A24): Evidence of Casualty Sensitivity Among Veterans and Non-veterans
#          Random Effects Specifications
########
m1 <- lmer(tone ~ mil*geo_fatal + rep + as.factor(date_month2) + (1| speakerid),  data=data)
m2 <- lmer(tone ~ mil*geo_fatal + rep + woman + age + 
             pip + obama + senate +
             word_count + word_count2  + as.factor(date_month2) +  (1| speakerid),  data = data)
m3 <- lmer(tone ~ mil*geo_fatal_cumulat + rep + as.factor(date_month2) + (1| speakerid),  data=data)
m4 <- lmer(tone ~ mil*geo_fatal_cumulat + rep + woman + age + 
             pip + obama + senate +
             word_count + word_count2  + as.factor(date_month2) +  (1| speakerid),  data = data)
stargazer(m1, m2, m3, m4,
          keep = c('mil','geo_fatal', 'geo_fatal_cumulat','nat_fatal_month','dwnom1','dwnom2','rep',
                   'woman', 'age', 'pip', 'obama', 'senate'), type="text")


########
# Table 4 (Table A25): Evidence of Casualty Sensitivity Among Veterans and Non-Vets, FE Specifications 
########
m1 <- lm(tone ~ geo_fatal +  word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0))
m2 <- lm(tone ~ geo_fatal  + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==0))
m3 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0))
m4 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==0))
m5 <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1))
m6 <- lm(tone ~ geo_fatal + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==1))
m7 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1))
m8 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==1))
stargazer(m1,m2, m3, m4, m5, m6, m7, m8,
          keep = c('geo_fatal','geo_fatal_cumulat','word_count','word_count2','months'),
          type="text")

stargazer(m1,m2, m3, m4, m5, m6, m7, m8,
          keep = c('geo_fatal','geo_fatal_cumulat','word_count','word_count2','months'))

#Quantities referenced in-text 
output1 <- summary(m1)
output2 <- summary(m3)
round(output1$coefficients[2]*max(data$geo_fatal),3) # -0.501, Max effect monthly
round(output2$coefficients[2]*max(data$geo_fatal_cumulat),3) # -0.379 Max effect cumulat
#Size of partisan gap
round(mean(data$tone[data$year==2006 & data$rep==1], na.rm=T),3)
round(mean(data$tone[data$year==2006 & data$rep==0], na.rm=T),3)
round(mean(data$tone[data$year==2006 & data$rep==1], na.rm=T) - mean(data$tone[data$year==2006 & data$rep==0], na.rm=T),3)

########
# Table 5: Sample Speeches and Speaker Tone
########
#Sample speeches and generate densities 
kennedy_sp <- subset(data, speech_id==1100027998) #-0.333
lewis_sp <- subset(data, speech_id==1080143807) #-3.2
graham_sp <- subset(data, speech_id==1130068681) #0.3
pdf("plots/tab5a.pdf",height=5,width=10)
par(mfrow=c(1,1), mar=c(6,2,13,2))
plot(NULL,
     xlim = c(-4,4), 
     ylim=c(0,0.3),
     axes = F, xlab = NA, ylab = NA) 
polygon(density(data$tone,bw=0.5)$x,
        density(data$tone,bw=0.5)$y ,
        col=alpha("steelblue",1),border="white",lwd=1.5)  
xlab("")
axis(side=1, las=1,cex.axis=2.5,mgp=c(3,1.5,0))
mtext("Tone of Speech", side = 1, line=3.75,cex=2.5) 
abline(v=kennedy_sp$tone, lwd=5, col="orange")
dev.off() 

pdf("plots/tab5b.pdf",height=5,width=10)
par(mfrow=c(1,1), mar=c(6,2,13,2))
plot(NULL,
     xlim = c(-4,4), 
     ylim=c(0,0.3),
     axes = F, xlab = NA, ylab = NA) 
polygon(density(data$tone,bw=0.5)$x,
        density(data$tone,bw=0.5)$y ,
        col=alpha("steelblue",1),border="white",lwd=1.5)  
xlab("")
axis(side=1, las=1,cex.axis=2.5,mgp=c(3,1.5,0))
mtext("Tone of Speech", side = 1, line=3.75,cex=2.5) 
abline(v=lewis_sp$tone, lwd=5, col="orange")
dev.off()

pdf("plots/tab5c.pdf",height=5,width=10)
par(mfrow=c(1,1), mar=c(6,2,13,2))
plot(NULL,
     xlim = c(-4,4), 
     ylim=c(0,0.3),
     axes = F, xlab = NA, ylab = NA) 
polygon(density(data$tone,bw=0.5)$x,
        density(data$tone,bw=0.5)$y ,
        col=alpha("steelblue",1),border="white",lwd=1.5)  
xlab("")
axis(side=1, las=1,cex.axis=2.5,mgp=c(3,1.5,0))
mtext("Tone of Speech", side = 1, line=3.75,cex=2.5) 
abline(v=graham_sp$tone, lwd=5, col="orange")
dev.off()

########
# Table A29: Veteran Status and the Use of Technical Military Terminology
########
m1 <- lmer(tech_econ_prop ~ mil + rep + woman + age   + pip + obama + senate + 
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),  data = data)
m2 <- lmer(tech_dod_prop ~ mil + rep + woman + age   + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),  data = data)
cov_list = c("Veteran", 
             "Republican", "Woman", 
             "Age", "President's Party", "Obama", 
             "Senate")

stargazer(m1, m2, 
          keep = c('mil', 'rep', 'woman',
                   'age',  'pip',  'obama', 'senate',
                   'word_count', 'word_count2'),
          covariate.labels = cov_list, type = 'text')

########
# Table A30: Cumulative Local Fatalitites and Uncertainty When Referencing Iraq
#          or Afghanistan
########
m1 <- lm(uncertainty1 ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0))
m2 <- lm(uncertainty1 ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==0))
m3 <- lm(uncertainty1 ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1))
m4 <- lm(uncertainty1 ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==1))
stargazer(m1, m2, m3, m4,  
          keep = c('geo_fatal_cumulat','word_count','word_count2','months'),
          type="text")

########
# Table A1: Veteran Status and Semantic Tone, Modified Lexicoder Dictionary
########
#Generate modified dictionary tone indicator
data$tone2 <- (data$lexi_modi_pos - data$lexi_modi_neg) / data$keyword_count 

m1 <- lmer(tone2 ~ mil + (1 | speakerid),  data=data)
m2 <- lmer(tone2 ~ mil + rep + (1 | speakerid),  data=data)
m3 <- lmer(tone2 ~ mil + rep + woman + age + 
             pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),  data = data)
m4 <- lmer(tone2 ~ mil + rep + mil_rep + woman + age + 
             pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),  data = data)
stargazer(m1, m2, m3, m4,
          keep = c('mil', 'rep', 'mil_rep', 'dwnom1', 'dwnom2','dw_mil', 'woman',
                   'age', 'geo_fatal_cumulat', 'pip', 'obama', 'senate',
                   'word_count', 'word_count2'), type="text")

stargazer(m1, m2, m3, m4,
          keep = c('mil', 'rep', 'mil_rep', 'dwnom1', 'dwnom2','dw_mil', 'woman',
                   'age', 'geo_fatal_cumulat', 'pip', 'obama', 'senate',
                   'word_count', 'word_count2'))


########
# Table A2: Evidence of Casualty Sensitivity Among Veterans and Non-Veterans, Random Effects Specifications, Modified Lexicoder Dictionary
########
m1 <- lmer(tone2 ~ mil*geo_fatal + rep + as.factor(date_month2) + (1| speakerid),  data=data)
m2 <- lmer(tone2 ~ mil*geo_fatal + rep + woman + age + 
             pip + obama + senate +
             word_count + word_count2  + as.factor(date_month2) +  (1| speakerid),  data = data)
m3 <- lmer(tone2 ~ mil*geo_fatal_cumulat + rep + as.factor(date_month2) + (1| speakerid),  data=data)
m4 <- lmer(tone2 ~ mil*geo_fatal_cumulat + rep + woman + age + 
             pip + obama + senate +
             word_count + word_count2  + as.factor(date_month2) +  (1| speakerid),  data = data)
stargazer(m1, m2, m3, m4,
          keep = c('mil','geo_fatal', 'geo_fatal_cumulat','nat_fatal_month','dwnom1','dwnom2','rep',
                   'woman', 'age', 'pip', 'obama', 'senate'), type="text")


########
# Table A3: Evidence of Casualty Sensitivity Among Veterans and Non-Veterans, Fixed Effects Specifications, Modified Lexicoder Dictionary
########
m1 <- lm(tone2 ~ geo_fatal +  word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0))
m2 <- lm(tone2 ~ geo_fatal  + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==0))
m3 <- lm(tone2 ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0))
m4 <- lm(tone2 ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==0))
m5 <- lm(tone2 ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1))
m6 <- lm(tone2 ~ geo_fatal + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==1))
m7 <- lm(tone2 ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1))
m8 <- lm(tone2 ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil==1))
stargazer(m1,m2, m3, m4, m5, m6, m7, m8,
          keep = c('geo_fatal','geo_fatal_cumulat','word_count','word_count2','months'), type='text')

########
# Table A22: Veteran Status and Semantic Tone (Combat Experience)
########
data$combat_rep <- data$combat * data$rep
data$mil_nocombat_rep <- data$mil_nocombat * data$rep
m1 <- lmer(tone ~ mil_nocombat + combat + (1 | speakerid),  data=data)
m2 <- lmer(tone ~ mil_nocombat + combat + rep + (1 | speakerid),  data=data)
m3 <- lmer(tone ~ mil_nocombat + combat + rep + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month)  + (1 | speakerid),  data = data)
m4 <- lmer(tone ~ mil_nocombat + mil_nocombat_rep + combat + combat_rep + rep + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month)  + (1 | speakerid),  data = data)
stargazer(m1, m2, m3, m4,
          keep = c('mil_nocombat', 'mil_nocombat_rep', 'combat', 'combat_rep', 'rep', 'woman', 'age', 
                   'geo_fatal_cumulat', 'pip', 'obama', 'senate', 'word_count', 'word_count2'), type="text")

stargazer(m1, m2, m3, m4,
          keep = c('mil_nocombat', 'mil_nocombat_rep', 'combat', 'combat_rep', 'rep', 'woman', 'age', 
                   'geo_fatal_cumulat', 'pip', 'obama', 'senate', 'word_count', 'word_count2'))
########
# Table A23: Monthly and Cumulative Local Fatalities and Semantic Tone When  
#           Referencing Iraq or Afghanistan by Veterans’ Combat Experience
########
m1 <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil_nocombat==1))
m2 <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,combat==1))
m3 <- lm(tone ~ geo_fatal + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil_nocombat==1))
m4 <- lm(tone ~ geo_fatal + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,combat==1))
m5 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil_nocombat==1))
m6 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,combat==1))
m7 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,mil_nocombat==1))
m8 <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + as.factor(date_month) + as.factor(speakerid),data=subset(data,combat==1))
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          type = "text", keep = c('geo_fatal', 'geo_fatal_cumulat', 
                                  'word_count','word_count2','months'))


########
# Table A26: National-Level Fatalities and Speaker Tone, Fixed Effects Specifications 
########
m1 <- lm(tone ~ nat_fatal_month +  word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0))
m2 <- lm(tone ~ log_nat_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0))
m3 <- lm(tone ~ nat_fatal_month + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1))
m4 <- lm(tone ~ log_nat_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1))
stargazer(m1, m2, m3, m4, type  = 'text',
          keep = c('nat_fatal_month', 'nat_fatal_cumulat', 
                   'log_nat_fatal_cumulat', 'log_nat_fatal',
                   'word_count','word_count2','months'))

########
# Table A27: National-Level Fatalities, Random Effects Specifications 
########
m1 <- lmer(tone ~ mil*nat_fatal_month + rep + months + (1| speakerid),  data=data)
m2 <- lmer(tone ~ mil*nat_fatal_month + rep + woman + age + pip +
             obama + senate +
             word_count + word_count2  + months +  (1| speakerid),  data = data)
m3 <- lmer(tone ~ mil*log_nat_fatal_cumulat + rep + months + (1| speakerid),  data=data)
m4 <- lmer(tone ~ mil*log_nat_fatal_cumulat + rep + woman + age + pip + 
             obama + senate +
             word_count + word_count2  + months +  (1| speakerid),  data = data)
stargazer(m1, m2, m3, m4,
          keep = c('mil','nat_fatal_month', 'nat_fatal_cumulat', 'log_nat_fatal_cumulat', 'log_nat_fatal', 'dwnom1','dwnom2','rep',
                   'woman', 'age', 'pip', 'obama', 'senate', 'months'), type="text")


########
# Table A28: Models with logged national cumulative casualties and local monthly casualties
########
m1 <- lmer(tone ~ mil + log_nat_fatal_cumulat + (1 | speakerid),  data=data)
m2 <- lmer(tone ~ mil + rep + log_nat_fatal_cumulat + (1 | speakerid),  data=data)
m3 <- lmer(tone ~ mil + rep + woman + age + 
             pip + obama + senate + log_nat_fatal_cumulat + 
             word_count + word_count2 + months  + (1 | speakerid),  data = data)
m4 <- lmer(tone ~ mil + rep + mil_rep + woman + age + 
             pip + obama + senate + log_nat_fatal_cumulat + 
             word_count + word_count2 + months  + (1 | speakerid),  data = data)
m5 <- lmer(tone ~ mil*geo_fatal + rep + months + log_nat_fatal_cumulat + word_count + word_count2  + (1| speakerid),  data=data)
m6 <- lmer(tone ~ mil*geo_fatal + rep + woman + age + pip + obama + senate +
             word_count + word_count2  + months + log_nat_fatal_cumulat +  (1| speakerid),  data = data)
m7 <- lm(tone ~ geo_fatal +  word_count + word_count2 + months + log_nat_fatal_cumulat + as.factor(speakerid),data=subset(data,mil==0))
m8 <- lm(tone ~ geo_fatal + word_count + word_count2 + months + log_nat_fatal_cumulat + as.factor(speakerid),data=subset(data,mil==1))
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          keep = c('mil','geo_fatal', 'geo_fatal',  'dwnom1','dwnom2','rep',
                   'woman', 'age', 'pip', 'obama', 'senate', 
                   'log_nat_fatal_cumulat' ,'months'),
          type="text")

########
# Table A31: Cumulative Local Fatalities and Uncertainty, Random Effects specifications
########
m1 <- lmer(uncertainty1 ~ mil*geo_fatal_cumulat + rep + as.factor(date_month2) + (1| speakerid),  data=data)
m2 <- lmer(uncertainty1 ~ mil*geo_fatal_cumulat + rep + woman + age + pip + 
             obama + senate +
             word_count + word_count2  + as.factor(date_month2) +  (1| speakerid),  data = data)
stargazer(m1, m2, 
          keep = c('mil','geo_fatal', 'geo_fatal',  'dwnom1','dwnom2','rep',
                   'woman', 'age', 'pip', 'obama', 'senate', 'months'), 
          type="text")

########
# Table A33: Veteran Status and Semantic Tone When Referencing 
#             the Abu Ghraib Scandal and the Troop Surge Strategy
########
# These spreadsheets just tag speeches in the relevant categories
ghraib <- read.csv( "speech_data/keyword_tags/ghraib_speech.csv")
surge <- read.csv("speech_data/keyword_tags/surge_speech.csv")
data <- merge(data, ghraib, by = "speech_id", all.x = T, all.y =F)
data$ghraib[is.na(data$ghraib)] <- 0
data <- merge(data, surge, by = "speech_id", all.x = T, all.y =F)
data$surge[is.na(data$surge)] <- 0
m1 <- lmer(tone ~  mil + ghraib  + dwnom1 + dwnom2 + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid) ,data=data)
m2 <- lmer(tone ~  ghraib*mil +mil + ghraib  + dwnom1 + dwnom2 + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid) ,data=data)
m3 <- lmer(tone ~ mil + surge + dwnom1 + dwnom2 + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),data=data)
m4 <- lmer(tone ~ surge*mil +mil + surge + dwnom1 + dwnom2 + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid),data=data)
stargazer(m1, m2, m3, m4, omit = "date_month2", type="text")

########
# Table A34: Veteran Status and Semantic Tone When Referencing 
#             Guantanamo Prison
########
gitmo2 <- read.csv("speech_data/keyword_tags/gitmo_speech2.csv")
data <- merge(data, gitmo2, by = "speech_id", all.x = T, all.y =F) 
data$gitmo2[is.na(data$gitmo2)] <- 0
m1 <- lmer(tone ~  mil + gitmo2  + dwnom1 + dwnom2 + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid) ,data=data)
m2 <- lmer(tone ~  gitmo2*mil +mil + gitmo2  + dwnom1 + dwnom2 + woman + age + 
             geo_fatal_cumulat + pip + obama + senate +
             word_count + word_count2 + as.factor(date_month2)  + (1 | speakerid) ,data=data)
stargazer(m1, m2, omit = 'date_month2', type="text")

########
# Table A35:Monthly Fatalities and Tone by Partisanship
########
crm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & rep==1))
cdm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & dem==1))
mrm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & rep==1))
mdm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & dem==1))
stargazer(crm, cdm,
          mrm, mdm,
          type="text",
          keep = c('geo_fatal','word_count','word_count2','months'))

########
# Table A36: Monthly Fatalities and Tone by Partisanship, Bush Admin.
########
crbm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & rep==1 & obama==0))
cdbm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & dem==1 & obama==0))
mrbm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & rep==1 & obama==0))
mdbm <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & dem==1 & obama==0))
stargazer(crbm, cdbm,
          mrbm, mdbm,
          keep = c('geo_fatal','word_count','word_count2','months'),
          type="text")

########
# Table A37: Monthly Fatalities and Tone by Partisanship, Obama Admin.
########
#OBAMA - Monthly Fatalities + Monthly Counter (A35)
crom <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & rep==1 & obama==1))
cdom <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & dem==1 & obama==1))
mrom <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & rep==1 & obama==1))
mdom <- lm(tone ~ geo_fatal + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & dem==1 & obama==1))
stargazer(crom, cdom,
          mrom, mdom,
          keep = c('geo_fatal','word_count','word_count2','months'), 
          type = "text")

########
# Table A38: Cumulative Fatalities and Tone by Partisanship, Full Sample
########
cr <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & rep==1))
cd <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & dem==1))
mr <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & rep==1))
md <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & dem==1))
stargazer(cr, cd,
          mr, md,
          keep = c('geo_fatal','word_count','word_count2','months'),
          type="text")

########
# Table A39: Cumulative Fatalities and Tone by Partisanship, Bush Admin
########
crb <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & rep==1 & obama==0))
cdb <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & dem==1 & obama==0))
mrb <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & rep==1 & obama==0))
mdb <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & dem==1 & obama==0))
stargazer(crb, cdb,
          mrb, mdb,
          keep = c('geo_fatal','word_count','word_count2','months'),
          type="text")

########
# Table A40: Cumulative Fatalities and Tone by Partisanship, Obama Admin
########
cro <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & rep==1 & obama==1))
cdo <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==0 & rep==0 & obama==1))
mro <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & rep==1 & obama==1))
mdo <- lm(tone ~ geo_fatal_cumulat + word_count + word_count2 + months + as.factor(speakerid),data=subset(data,mil==1 & dem==1 & obama==1))
stargazer(cro, cdo,
          mro, mdo,
          keep = c('geo_fatal','word_count','word_count2','months'),
          type="text")

########
# Figure A4: Casualty Sensitivity Among Partisan Subgroups, Monthly Fatalities
########
pdf("plots/figA4.pdf",height=6,width=12)
par(mar=c(4,4,2,2),
    mfrow=c(1,3))
coefs <- c(coefficients(crm)[2], coefficients(cdm)[2],
           coefficients(mrm)[2], coefficients(mdm)[2])
ses <- c(sqrt(vcov(crm)[2,2]), sqrt(vcov(cdm)[2,2]), 
         sqrt(vcov(mrm)[2,2]), sqrt(vcov(mdm)[2,2]))
plot(NULL,
     ylim = c(1,length(coefs)+.05), 
     xlim=c(-.14,0.14),
     axes = F, xlab = NA, ylab = NA) 
abline(v=0, col="#b2182b",lty=5, lwd=1)
for(ii in 1:length(coefs)){
  lines(c(coefs[ii]-(1.96*ses[ii]),coefs[ii]+(1.96*ses[ii])),
        c(ii,ii),
        col="black", lwd=1.75)
}
points(coefs,seq(1,length(coefs),1),bg=alpha("gray80",1), col=alpha("black",1),lwd=1.75,cex=1, pch=21)
text(coefs,seq(1,length(coefs),1),
     labels=c("Non-Veteran, Republican", "Non-Veteran, Democrat", "Veteran, Republican", "Veteran, Democrat"),
     cex=1.2, pos=3)
axis(1)
mtext(side = 1, "Coefficient Value", line =2.5,cex=1) # label bottom axis
title("Full Sample", line=1)
box()

coefs <- c(coefficients(crbm)[2], coefficients(cdbm)[2],
           coefficients(mrbm)[2], coefficients(mdbm)[2])
ses <- c(sqrt(vcov(crbm)[2,2]), sqrt(vcov(cdbm)[2,2]), 
         sqrt(vcov(mrbm)[2,2]), sqrt(vcov(mdbm)[2,2]))
plot(NULL,
     ylim = c(1,length(coefs)+.05), 
     xlim=c(-.14,0.14),
     axes = F, xlab = NA, ylab = NA) 
abline(v=0, col="#b2182b",lty=5, lwd=1)
for(ii in 1:length(coefs)){
  lines(c(coefs[ii]-(1.96*ses[ii]),coefs[ii]+(1.96*ses[ii])),
        c(ii,ii),
        col="black", lwd=1.75)
}
points(coefs,seq(1,length(coefs),1),bg=alpha("gray80",1), col=alpha("black",1),lwd=1.75,cex=1, pch=21)
text(coefs,seq(1,length(coefs),1),
     labels=c("Non-Veteran, Republican", "Non-Veteran, Democrat", "Veteran, Republican", "Veteran, Democrat"),
     cex=1.2, pos=3)
axis(1)
mtext(side = 1, "Coefficient Value", line =2.5,cex=1) # label bottom axis
title("Bush Administration", line=1)
box()
coefs <- c(coefficients(crom)[2], coefficients(cdom)[2],
           coefficients(mrom)[2], coefficients(mdom)[2])
ses <- c(sqrt(vcov(crom)[2,2]), sqrt(vcov(cdom)[2,2]), 
         sqrt(vcov(mrom)[2,2]), sqrt(vcov(mdom)[2,2]))
plot(NULL,
     ylim = c(1,length(coefs)+.05), 
     xlim=c(-.3,.3),
     axes = F, xlab = NA, ylab = NA) 
abline(v=0, col="#b2182b",lty=5, lwd=1)
for(ii in 1:length(coefs)){
  lines(c(coefs[ii]-(1.96*ses[ii]),coefs[ii]+(1.96*ses[ii])),
        c(ii,ii),
        col="black", lwd=1.75)
}
points(coefs,seq(1,length(coefs),1),bg=alpha("gray80",1), col=alpha("black",1),lwd=1.75,cex=1, pch=21)
text(coefs,seq(1,length(coefs),1),
     labels=c("Non-Veteran, Republican", "Non-Veteran, Democrat", "Veteran, Republican", "Veteran, Democrat"),
     cex=1.2, pos=3)
axis(1)
mtext(side = 1, "Coefficient Value", line =2.5,cex=1) # label bottom axis
title("Obama Administration", line=1)
box()
dev.off()


########
# Figure A5: Casualty Sensitivity Among Partisan Subgroups, Cumulative Fatalities
########
pdf("plots/figA5.pdf",height=6,width=12)
par(mar=c(4,4,2,2),
    mfrow=c(1,3))
coefs <- c(coefficients(cr)[2], coefficients(cd)[2],
           coefficients(mr)[2], coefficients(md)[2])
ses <- c(sqrt(vcov(cr)[2,2]), sqrt(vcov(cd)[2,2]), 
         sqrt(vcov(mr)[2,2]), sqrt(vcov(md)[2,2]))
plot(NULL,
     ylim = c(1,length(coefs)+.05), 
     xlim=c(-.005,.005),
     axes = F, xlab = NA, ylab = NA) 
abline(v=0, col="#b2182b",lty=5, lwd=1)
for(ii in 1:length(coefs)){
  lines(c(coefs[ii]-(1.96*ses[ii]),coefs[ii]+(1.96*ses[ii])),
        c(ii,ii),
        col="black", lwd=1.75)
}
points(coefs,seq(1,length(coefs),1),bg=alpha("gray80",1), col=alpha("black",1),lwd=1.75,cex=1, pch=21)
text(coefs,seq(1,length(coefs),1),
     labels=c("Non-Veteran, Republican", "Non-Veteran, Democrat", "Veteran, Republican", "Veteran, Democrat"),
     cex=1.2, pos=3)
axis(1)
mtext(side = 1, "Coefficient Value", line =2.5,cex=1) # label bottom axis
title("Full Sample", line=1)
box()
coefs <- c(coefficients(crb)[2], coefficients(cdb)[2],
           coefficients(mrb)[2], coefficients(mdb)[2])
ses <- c(sqrt(vcov(crb)[2,2]), sqrt(vcov(cdb)[2,2]), 
         sqrt(vcov(mrb)[2,2]), sqrt(vcov(mdb)[2,2]))
plot(NULL,
     ylim = c(1,length(coefs)+.05), 
     xlim=c(-.005,0.01),
     axes = F, xlab = NA, ylab = NA) 
abline(v=0, col="#b2182b",lty=5, lwd=1)
for(ii in 1:length(coefs)){
  lines(c(coefs[ii]-(1.96*ses[ii]),coefs[ii]+(1.96*ses[ii])),
        c(ii,ii),
        col="black", lwd=1.75)
}
points(coefs,seq(1,length(coefs),1),bg=alpha("gray80",1), col=alpha("black",1),lwd=1.75,cex=1, pch=21)
text(coefs,seq(1,length(coefs),1),
     labels=c("Non-Veteran, Republican", "Non-Veteran, Democrat", "Veteran, Republican", "Veteran, Democrat"),
     cex=1.2, pos=3)
axis(1)
mtext(side = 1, "Coefficient Value", line =2.5,cex=1) # label bottom axis
title("Bush Administration", line=1)
box()
coefs <- c(coefficients(cro)[2], coefficients(cdo)[2],
           coefficients(mro)[2], coefficients(mdo)[2])
ses <- c(sqrt(vcov(cro)[2,2]), sqrt(vcov(cdo)[2,2]), 
         sqrt(vcov(mro)[2,2]), sqrt(vcov(mdo)[2,2]))
plot(NULL,
     ylim = c(1,length(coefs)+.05), 
     xlim=c(-.05,0.025),
     axes = F, xlab = NA, ylab = NA) 
abline(v=0, col="#b2182b",lty=5, lwd=1)
for(ii in 1:length(coefs)){
  lines(c(coefs[ii]-(1.96*ses[ii]),coefs[ii]+(1.96*ses[ii])),
        c(ii,ii),
        col="black", lwd=1.75)
}
points(coefs,seq(1,length(coefs),1),bg=alpha("gray80",1), col=alpha("black",1),lwd=1.75,cex=1, pch=21)
text(coefs,seq(1,length(coefs),1),
     labels=c("Non-Veteran, Republican", "Non-Veteran, Democrat", "Veteran, Republican", "Veteran, Democrat"),
     cex=1.2, pos=3)
axis(1)
mtext(side = 1, "Coefficient Value", line =2.5,cex=1) # label bottom axis
title("Obama Administration",line=1)
box()
dev.off()













